home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1995 April / macformat-023.iso / Shareware City / Developers / RLaB / rlib / compan.r < prev    next >
Encoding:
Text File  |  1994-09-23  |  1.8 KB  |  70 lines  |  [TEXT/RLAB]
  1. //-------------------------------------------------------------------//
  2.  
  3. // Synopsis:    Companion matrix.
  4.  
  5. // Syntax:    compan ( P )
  6.  
  7. // Description:
  8.  
  9. //    There are three cases.
  10.  
  11. //    If P is a scalar then COMPAN(P) is the P-by-P matrix COMPAN(1:P+1).
  12. //      If P is an (n+1)-vector, COMPAN(P) is the n-by-n companion matrix
  13. //           whose first row is -P(2:n+1)/P(1).
  14. //      If P is a square matrix, COMPAN(P) is the companion matrix
  15. //           of the characteristic polynomial of P, computed as
  16. //           COMPAN(POLY(P)).
  17.  
  18. //      References:
  19. //      J.H. Wilkinson, The Algebraic Eigenvalue Problem,
  20. //           Oxford University Press, 1965, p. 12.
  21. //      G.H. Golub and C.F. Van Loan, Matrix Computations, second edition,
  22. //           Johns Hopkins University Press, Baltimore, Maryland, 1989,
  23. //           sec 7.4.6.
  24. //      C. Kenney and A.J. Laub, Controllability and stability radii for
  25. //          companion form systems, Math. Control Signals Systems, 1 (1988),
  26. //          pp. 239-256. (Gives explicit formulas for the singular values of
  27. //          COMPAN(P).)
  28.  
  29. //    This file is a translation of compan.m from version 2.0 of
  30. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  31. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  32.  
  33. //-------------------------------------------------------------------//
  34.  
  35. rfile poly
  36.  
  37. compan = function ( p )
  38. {
  39.   local (p)
  40.  
  41.   m = p.nr; n = p.nc;
  42.  
  43.   if (n == m && n > 1)
  44.   {
  45.     // Matrix argument.
  46.     return $self (poly (p));
  47.   }
  48.  
  49.   n = max (n,m);
  50.   //  Handle scalar p.
  51.   if (n == 1)
  52.   {
  53.     n = p+1;
  54.     p = 1:n;
  55.   }
  56.  
  57.   p = p[:]';                    // Ensure p is a row vector.
  58.  
  59.   // Construct matrix of order n-1.
  60.   if (n == 2)
  61.   {
  62.     A = 1;
  63.   else
  64.     A = diag(ones(1,n-2),-1);
  65.     A[1;] = -p[2:n]/p[1];
  66.   }
  67.  
  68.   return A;
  69. };
  70.